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Pattern formation exhibited by a two-dimensional reaction-diffusion system 



^ ■ in the fast inhibitor hmit is considered from the point of view of interface 

^^ ■ motion. A dissipative nonlocal equation of motion for the boundary between 

f~^ ■ high and low concentrations of the slow species is derived heuristically. Under 

^: 

0^ . these dynamics, a compact domain of high concentration may develop into a 



space-filling labyrinthine structure in which nearby fronts repel. Similar pat- 
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i-rt ■ terns have been observed recently by Lee, McCormick, Ouyang, and Swinney 

O ■ in a reacting chemical system. 
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In the study of chemical systems with both reactions and diffusion, one may discern two 
broad classes of spatial patterns: extended and compact. Extended patterns typically arise 
from supercritical symmetry-breaking bifurcations |1|. In the two dimensional case, which 
we consider here, they are often regular, periodic structures such as arrays of stripes, discs, 
or hexagons 0. Compact patterns, or localized states, appear in systems with subcritical 
bifurcations via a nucleation process ^, and typically take the form of a single one of the 
repeating units found in extended systems. Both classes of patterns are often described in 
terms of a competition between two chemical species: an autocatalytic "activator" and its 
"inhibitor." 

Localized states can exhibit a fingering instability [^ . One can imagine that these fingers 
may grow and branch until a complicated labyrinthine pattern fills the entire plane. Such 
structures may actually be metastable states; barriers to domain fission may prevent the 
pattern from evolving into the ground state (presumably a regular array of stripes) from 
an initial condition which is topologically different. Qualitatively similar kinds of pattern 
formation appear in other systems |^. 

Complicated, labyrinthine pattern evolution in a chemical system has recently been ob- 
served by Lee, McCormick, Ouyang, and Swinney ^ in an iodate-ferrocyanide-sulfite reac- 
tion (see Fig. |l]). The patterns are composed of regions with one of two different chemical 
compositions. This system appears to be bistable; if it is prepared with one of the two possi- 
ble uniform compositions, it will persist in that state. Nontrivial pattern formation requires 
a nucleation site in an otherwise uniform background. These experiments also indicate that 
in the formation of patterns, nearby boundaries or fronts repel each other. We will use these 
results as a guide to the general features one would like to see in a model system 0. 

It is natural to seek a representation of the dynamics of these and similar patterns in 
terms of the interface between domains of different composition; in order to capture the 
properties of repulsion and nonintersection, such a contour dynamics must be nonlocal, 
coupling segments of the interface that are distant in arclength yet close in space. The 
purpose of this letter is twofold: to give an intuitive construction of a fully nonlocal curve 



dynamics from a set of react ion- diffusion equations and to show that it is useful for describing 
pattern formation. The evolution equation is written in terms of the intrinsic geometry of the 
curve, and is valid for pattern evolution far beyond the linear instability of localized states. 
In this simplified form, it is easier to identify the physics responsible for the destabilization of 
a compact pattern, the tendency of a pattern to grow or shrink, and the interaction between 
different portions of the interface. The curve dynamics is studied numerically and is shown 
to reproduce the quahtative features of the experimental patterns. 

The reaction-diffusion pair studied here is similar to models of spiral wave formation p| 
and nerve impulse propagation ^. With u the activator and v the inhibitor, we consider 

ut = DS/'^u - Y\u) - pv , (la) 

evt = V^w — av + Pu . (lb) 

The nonlinear function F(M),with derivative F', embodies the autocatalytic nature of the 
activator. It is typically a polynomial in u with a double well structure whose minima, not 
necessarily of equal depth, we label u±. Patterns like that in Fig. |I| can then form in which 
regions of m ^ m+ (e.g. white) are surrounded by a region of m ^ m_ (black). Inhibition 
of u is achieved for p > 0, while the couplings a and (3 reflect the self-limiting behavior of 
the inhibitor and its stimulation by the activator, respectively. The small number e defines 
the fast-inhibitor limit. This limit is opposite to the limit assumed in phase-field models 
||T0| and spiral wave dynamics Pj. The fast inhibitor assumption appears to be the simplest 
assumption that allows the elimination of one of the fields, giving rise to spatial nonlocality 
for the remaining field without introducing temporal nonlocality. A similar calculation can 



be done for the slow-inhibitor limit [jTT|, but due to temporal nonlocality, the resulting 
equations are more complicated, obscuring the essential physics. 

We begin by discussing the interaction of fronts, illustrating the property of self-avoidance 
present in the reaction-diffusion pair (1), for certain parameter ranges. Figure 2 shows a 
simulation of a one-dimensional version of (1) with periodic boundary conditions. The 



patterns are defined by the sharp interfaces of the u field. The parameters are such that 
the u ^ n_|_ region expands into the u c::^ u_ region. When the interfaces get too close, the 
exponential tails of the v field begin to overlap, causing a repulsion which stabilizes the final 
field configuration. In the derivation of the two-dimensional contour dynamics, we shall 
analytically see the source of this repulsion. 

To derive an interface evolution equation for two-dimensional patterns, we take the fast- 
inhibitor limit, setting e = and thus slaving v to u, 

t;(x,t) = /3y'dV^(|x-x'|)M(x') (2) 

where the Green's function G{r) = Ko(r/^)/(27r), and ^ = a~^/^. Substituting for v in 



equations (pla]), we obtain the variational dynamics Ut = —5J-'[u]/6u , with 

+ pI fd^xfd^x'u{^)g{\^ - x'l)n(x'). (3) 

The functional JF decreases monotonically in time, possibly reaching a local minimum. For 
e > 0, there is no such JF. The necessity of a large- amplitude perturbation for the initiation of 
nontrivial pattern formation follows directly from the bistability of F{u) and the variational 
form of the dynamics. 

We would like to find an equation of motion for the interface r(s) in the form rj = 
—6J-'[r]/6r, starting from the equation Ut = —6J-'[u]/6u . To determine the boundary func- 
tional JF[r] from J-'[u] in (|^), we restrict ourselves to studying the dynamics of a single island 
of Uin — M+ in an infinite sea of Uout — u^. We make two assumptions to evaluate the term 
involving Vm: the profile is sharp, so that the gradients are localized on the boundary F, 
and the curvature of the boundary does not significantly affect the interface profile. The 
derivation of the contour dynamics using asymptotic expansions will be discussed elsewhere 



11|. We proceed here with a heuristic derivation because it yields greater physical insight. 



The integral of (Vu) is then proportional to the perimeter of the shape, the constant of 
proportionality being approximately D{uin — Uout)^/l, with / a characteristic length scale of 
the profile. One cannot determine / further without specifying the precise form of F. 



To evaluate the remaining terms in JF, we take u to be piecewise constant. The first 
term gives jF{u) = F{uin) /j„ + F{uout) lout ' where in and out refer to the areas inside 
and outside of F. The term involving Q can be written as 



ffngu= ulf f g + ul, f f g 

J J J inJ in JoutJout 



+2UinUout g ■ (4) 

JinJout 

This may be further simplified by an application of Stokes' theorem and the defining relation 
for g, (V^ — ^^'^) ^(x) = — (5(x). Apart from an unimportant constant, the energy is then 

J^[r] = UA + -fL- f^ ids jds'l ■ i'g (R) , (5) 

where 7 ~ DA'^/21, 

n ^ F(«,„) - F{uout) + {pPe/2) {ul - «L) , (6) 

and A = Um — Uout, R = r(s) — r(s'), R = |R|. We interpret 7 as a line tension (associated 
with the boundary length L) and 11 as a pressure (associated with the enclosed area A). 
The pressure, which can be of either sign, reflects the difference between the values of 
the energy F(u) inside and outside F. Note from the definition of FI that the coupling 
of u to f has created an effective energy F{u) + {p/3C,'^/2)u'^ whose minima determine the 
values Uin and Uout- The nonlocal terms in (^ (and (0), below) have the appearance of 
a self-induction interaction, although g provides an exponential screening. While at first 
sight somewhat unusual, a coupling between tangent vectors appears in other systems with 
piecewise constant fields |1T2|J13[| . 

Next we determine the equation of motion of r from that of u. This has two parts: (i) 
relating the time derivatives Ut and r^, and (ii) relating the functional derivatives 6/6u and 
6/6r. Since Ut is nonzero only near the boundary, the approximation Ut — {A/l)h-rt is valid. 
Likewise, the variation 6/6u is large only near the boundary, suggesting the identification 
jTSl {6/6u) -^ (//A)n ■ {S/6r) . Rescaling time by (A//)~ , the equation of motion is 

h ■ Tt = — n — 7k(s) 



-p(3C^ fds'R{s, s') X t(s') g'{R). (7) 

The prime on Q indicates a derivative with respect to R, and the cross product is a scalar 
in two dimensions. 

Note that the repulsion between adjacent fronts depends only on the fact that Q > 0, 
and not on the specific form of Q (as can be seen from (|^) and the variational form of the 
equation of motion). This is important, because in a more exact derivation of the contour 
dynamics that takes the interface profile into account, the function appearing in equation 
(^) may not be Q, but rather one with a less singular behavior as -R — * fl^. Note however 
that Q has only a logarithmic, integrable singularity at the origin, so the dynamics derived 
here is well defined. Even in the presence of such a cutoff in Q the dynamics would not 
be qualitatively affected, since the sign of Q alone determines whether fronts are attracted 
or repelled and the large R behavior of the purported function would be unchanged (i.e., 
exponential decay). 

By further redefining the time scale such that the coefficient of the nonlocal term in (|^ 
is unity, one finds only three relevant parameters remain: ^ and a rescaled 7 and 11. The 
model (0) was investiged numerically with a pseudospectral technique. Figure ^ shows the 
evolution of a circle seeded with sinusoidal perturbations. A complicated labyrinth forms 
due to the repeated fingering of the boundary. The interface never crosses itself, due to the 
repulsion between the adjacent fingers. The competition between the inward pressure 11 and 
the repulsion within a finger sets the finger width, while the interfinger distance is set by 
the inhibitor length scale ^ alone. The nonlocal nature of the dynamics severely limits the 
time scale over which the evolution may be followed numerically; with further computation, 
the pattern in Fig. ^ would continue to evolve beyond the final picture. 

Figure ^ shows the evolution of an island with a larger value of 11 than that in Fig. |^, 
implying a larger energy difference between the two possible homogeneous states of u. The 
area enclosed by the interface changes dramatically over time, the pattern simply shrinking 
to a circle. This shrinkage of a fingered structure has also been seen in the work of Lee, et 
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a/, i. 

To gain insight into the basic mechanism of the fingering instabihty, we develop an 
approximate local dynamics valid when ^k ^ 1 and the shape is approximately circular 
p!5| . Observe that the most important contribution to the double integral in (|^) is from 
the region |s' — s| < ^. Near the point r(s), expand the scalar product as t(s) • t(s') ~ 
1 — (l/2)(s — s')^/t(s)^ + ■ ■ ■ 5 then perform the s' integral over the extended region [— oo, oo] 
to obtain an effective local energy functional JF„: 

^„ ~ j^ds |7. + l^i^^Y + O [{iKf) I . (8) 

The first term contributes to an effective line tension, while the second we recognize from 
elasticity theory as the bending energy of a rigid rod. A detailed calculation yields 7^, = 
— p^^/4 and k^ = p^^/8. While k^ is positive, 7^ is negative. The effective line tension 
7 = 7 + 7^, can be negative, favoring proliferation of the interface. Again, note that the sign 
of 7^ is opposite that of Q, and hence the physics of fingering depends only on the sign of Q 
and not its specific functional form. 

Under the local approximation, the interface motion is 

-n - 7K - fc,e' f/tss + 2'^') ' (9) 



similar to "geometric" models |]T6| of crystal growth. In a linear stability analysis of a 
circular shape, the growth rate cr„ of the n^^ mode is cr„ ~ —777,^ — k^C.'^n'^. For 7 < 0, 
there is a band of unstable modes whose maximum extent is limited by the rigidity. In the 
numerical studies shown here, the effective tension is negative in both cases. In Fig. ^, the 
inward pressure is insufficient to prevent the proliferation of the interface, while in Fig. ^ 
the larger pressure overcomes the negative surface tension. 

In summary, starting from a reaction-diffusion system in two dimensions, we have con- 
structed a boundary dynamics useful in describing pattern formation. Possible extensions 
include studying inertial effects from higher order terms in e, and the application of these 
methods to higher-dimensional problems. It is also of interest to determine if the partic- 
ular chemical kinetics in the experiments of Lee, et al. may indeed be mapped onto the 



interface model developed here. Finally, it is possible to extend these methods to study 
the more physically relevant slow-inhibitor problem |Tl|. Perhaps issues such as spiral wave 
stability @|T^ could be given a more intuitive interpretation when studied from a geometric 
viewpoint. 

We are grateful to K.J. Lee, W.D. McCormick, Q. Ouyang, and H.L. Swinney for com- 
munication of their results prior to publication and for kindly providing Fig. |l|, and thank 
V. Hakim and S. Leibler for important discussions on localized states, and D. Barkley, S. 
Erramilli, M.J. Shelley, and A.T. Winfree for discussions and correspondence. This ma- 
terial is based upon work supported under an N.S.F. Graduate Fellowship (DMP), N.S.F. 
Presidential Faculty Fellowship Grant DMR-9350227 and the Alfred P. Sloan Foundation 
(REG). 
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FIGURES 
FIG. 1. Pattern formation in the chemical system of Lee, McCormick, Ouyang, and Swinney 

Ig]. White and black regions correspond, respectively, to low and high pH, made visible with a pH 

indicator. Figure courtesy of Lee, et al. [6]. 

FIG. 2. Space-time portrait of the interaction of two fronts, from numerical solution of Eq. 1 
with e = 0.008,L> = 0.2, p = 0.1, a = (3 = 0.2, F'{u) = au + hu^ + dn^, with a = 0.23, h = -1.23, 
and c = 1.0. Solid lines show M(x,t), dashed lines are w(x,t), with time increasing upward. 

FIG. 3. Numerical solution of (j^) with 7 = 2.0, 11 = 0.02, and ^ = 1.0. The initial condition is 
a perturbed circle of radius 15, and the time interval between shapes is 12. 

FIG. 4. Shrinkage of a domain with 7 = 2.0, 11 = 2.0, and ^ = 1.0. The time interval between 
shapes is 1. 
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